Statistical mechanics model of angiogenic tumor growth 
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We examine a lattice model of tumor growth where survival of tumor cells depends on the sup- 
plied nutrients. When such a supply is random, the extinction of tumors belongs to the directed 
percolation universality class. However, when the supply is correlated with distribution of tumor 
cells, which as we suggest might mimick the angiogenic growth, the extinction shows different, and 
most likely novel critical behaviour. Such a correlation affects also the morphology of the growing 
tumors and drastically raise tumor survival probability. 
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Due to series of mutations of genes responsible for pro- 
liferation, some cells might initiate the abnormal and un- 
controlled growth process, commonly named as tumor 
In the first stage of the process, called avascular 
growth, there is no blood supply and the growth is limited 
by the amount of oxygen and nutrients that the tumor 
receives through its surface. The growing tumor initiates, 
however, a number of accompanying processes in the host 
environment. In particular, large demand for nutrients 
stimulates tumor cells to produce angiogenic factors that 
regulate the formation and growth of new blood vessels 
in the region. This stage is called angiogenesis and it in- 
termediates between avascular and vascular growth. The 
vascular growth is the third stage, which begins when 
the blood vessels have reached the tumor. In this stage 
the tumor receives a vast amount of nutrients and can 
grow much larger than it was possible during the avas- 
cular growth. Moreover, the vasculature might be used 
to spread tumor cells throughout the body of the host, 
which very often leads to its death. 

Tumor growth is a very complex process and to fully 
understand its nature, one has to resort to computa- 
tional techniques, which would supplement biological and 
medical approaches Various models were used to 

describe tumor growth 7|. Initially, they were contin- 
uous models formulated in terms of partial differential 
equations and studied mainly from mathematical point 
of view 0. More recently, inclusion of biomechanical 
details and coupling of tumor growth with the vascular- 
ization process shifted modeling toward more physical 
approach '5]. 

An alternative to differential-equation approach is 
based on discrete lattice models such as, for example, 
nonequilibrium Q-state Potts models. With such a mod- 
eling, one can implement several aspects of cell dynamics, 
which are difficult to treat simultaneously using continu- 
ous modeling, as e.g. multiplication, competition, aging, 
death, mutations and even adhesion or chemotaxis [lOj. 
Some models of this kind were successfully tested against 
clinical data of certain forms of cancer 11, 131. One can 



hope that further development will result in lattice mod- 
els taking into account some other important, but so far 
neglected, factors like heterogeneity, immune cells or the 
role of chemoattractants. 

Discrete lattice modeling allows us to describe tumor 
growth problem using tools developed in statistical me- 
chanics for studying complex systems. A distinctive fea- 
ture of these systems is spontaneous emergence of cer- 
tain properties, which cannot be traced to the character 
of individual parts. It is believed that life, conscious- 
ness, or functioning of ant colonies are examples of such 
emergent phenomena [l3| . There is an increasing evi- 
dence that cancer can also be considered as an emergent 
property, and thus, developing statistical-mechanics ap- 
proaches seems to be very promising 14-3]. In partic- 
ular, one can hope that such bottom-up modeling will 
help us to explain tumor growth in terms of cell parame- 
ters, which might contribute to its better prediction and 
control. 

The ultimate goal, namely builiding realistic and 
testable against real data models, will most likely re- 
quire very complex, multi-scale models, which will be 
difficult to understand without extensive computer sim- 
ulations. To develop some intuitive understanding of tu- 
mor growth, it is thus desirable to examine some simpler 
models, which hopefully contain important ingredients 
of the process. In the present Letter, we examine a sim- 
ple lattice growth model where tumor cells survival and 
breeding depend on the supplied nutrients. Mimicking 
angiogenetic processes, we assume that the supply of nu- 
trients is positively correlated with the distribution of 
tumor cells. It turns out that such a correlation substan- 
tially changes the statistical mechanics behaviour of the 
model. In particular, the extinction of tumor does not 
belong to the expected directed-percolation universality 
class and the critical exponent (3 describing the order pa- 
rameter most likely gets the classical value /? = 1, even in 
the d = 1 version of the model. Such a correlation affects 
also the morphology of growing tumors and drastically 
raises the tumor-survival probability. 
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In our model, each site of a d-dimensional lattice either 
is occupied by tumor, nutrient, both tumor and nutri- 
ent, or is empty. At a rate p, nutrients are supplied to 
a chosen site of the lattice, provided that the site is not 
occu pied already by a nutrient. The roulette- wheel selec- 
tion [17[ is used to choose the site for such a supply and 
the corresponding weight w depends on whether the site 
is occupied by tumor (w = 1 -I- A) or not {w = 1). The 
parameter A > takes into account angiogenic effects 
of increased nutrients supply due to formation of new 
blood vessels in the vicinity of tumor cells. At a rate 
1 — p, a tumor cell on a randomly chosen site is updated. 
The tumor cell survives provided that there is a nutri- 
ent on this site, otherwise it dies. The surviving tumor 
cell consumes the nutrient and attempts to breed pro- 
vided that there is a site without a tumor cell among its 
nearest neighbours. Closely related models but without 
any preference in nutrient supply (A = 0) were already 
studied 

To examine our model we used Monte Carlo simula- 
tions [31 . For various p and A, we measured the steady- 
state densities of tumor cells Xt and of nutrients Xn- Sim- 
ulations started from a random initial configuration and 
the model relaxed until a steady state was reached. We 
examined lattices of various sizes N to ensure that the 
obtained results are TV- independent. We also measured 
the time dependence of the tumor cell density xt (t) , with 
the unit of time corresponding to N update attempts. 
The density xt{t) for each t is an average over indepen- 
dent runs. At the critical point, the density xt{t) is ex- 
pected to have a power-law decay Xt{t) ~ , where S is 
a characteristic exponent (20j . 

First, we describe results for the d = 1 version of the 
model. Simulations show that, for sufficiently large p, 
the model remains in an active phase with xt > 0, which 
terminates at a critical point Pc depending on A. For 
p < Pc, the steady state of the model is an absorb- 
ing state Xt — and Xn = 1 (tumor cells die out of 
the lack of nutrients). At p — pc, the model undergoes 
the phase transition from an active into absorbing phase 
and we expect that xt is the corresponding order pa- 
rameter. As it is already known, models with a single 
absorbing state are expected to belong to the so-called 
directed-percolation (DP) universality class [lOl- In this 
universality class and for d ^ 1, the critical exponent 
Sdp — 0.159(1) and the decay of the order parameter 
upon approaching the critical point is described by the 
exponent /3_dp — 0.276(1). The calculated steady-state 
values of Xt and Xn for A = and 3 are presented in 
FiglH For A = 0, we estimate Pc = 0.3691(2) and the 
exponent f3 is very close to the DP value (inset in FigllJ. 
The decay of Xt{t) at the critical point (Figl2]) confirms 
the DP universality class in this case. 

Much different behaviour is seen for A = 3. In this 
case, we estimate Pc = 0.1765(2) and the exponent 
/? = 1.0(1) (inset in FigH]). The same estimate of pc 
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FIG. 1. Steady-state densities of nutrients (xn) and of tumor 
cells {xt) as a function of p calculated for d = 1 and A = 
and 3. The inset shows the behaviour of Xt close to the critical 
point (log-log scale). While for A = the directed-percolation 
scaling is seen, much different behaviour appears for A = 3. 




FIG. 2. Time dependence of the density of tumor cells xt{t) 
(log-log scale) for d — 1, A = and (from top) p = 0.37, 
0.3695, 0.3692, 0.36905, 0.3689 0.368, 0.366, and 0.36. The 
dotted straight line has the slope corresponding to the di- 
rected percolation value 5 = 0.159. 



is obtained from the behaviour of xt{t) (Fig. |3]). How- 
ever, at p — Pc the decay seems to be described by the 
exponent 5 — 0.60(5), which is again much different than 
the DP value 0.159. In calculations for A = and 3, and 
close to critical points, systems as large as = 5x10*^ 
were used and simulation and relaxation times were of 
the order of 10^. Further away from critical points, less 
extensive simulations were required. 

For A = 0.1 on a shorter time scale {t 10^), one can 
notice a slower, DP-like decay of Xt{t), that, however, 
turns to much faster decay on a longer time scale (FigH)) . 
Such a behaviour indicates the proximity of the DP be- 
haviour that occurs at A = 0. However, critical expo- 
nents /? and S for A = 0.1 are nearly the same as for 
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FIG. 3. Time dependence of the density of tumor cells Xt{t) 
(log-log scale) for d = 1, A = 3 and (from top) p = 0.19, 0.18, 
0.178, 0.177, 0.1765, 0.176, 0.175, 0.17. The dotted straight 
line has the slope corresponding to S — 0.6. 



A = 3. We do not present numerical results, but simula- 
tions show that the same estimations of critical exponents 
are obtained for A = 5 and 10. Hence, our results suggest 
that the critical exponents /3 = 1-0(1) and S = 0.60(5) 
are universal most likely for any A > 0. 

Let us notice that for increasing A, the critical rate pc 
decreases and the tumor-free phase shrinks. Thus, an- 
giogenic factors, that in our model correspond to positive 
correlation between nutrient supply and tumor cell distri- 
bution (A > 0) make starvation of tumors to death more 
difficult. Such a behaviour of our model is very plau- 
sible since angiogenesis was "invented" by tumors just 
for this particular reason. For very large A, the critical 
rate Pc seems to vanish, which means that tumors can 
survive even under a very small rate of the nutrient sup- 
ply. Although the presented model is too simple to de- 
scribe detailed complexity of real tumors, in our opinion, 
it qualitatively correctly captures the role of angiogenic 
factors. 

Despite possessing a single absorbing state, our model 
most likely for any A > does not belong to the directed- 
percolation universality class. The exponent /3 seems 
to take the classical, mean-field value (3 — 1 that for 
a, d = 1 model is certainly a puzzling result. The 
exponent S = 0.60(5) is different than the mean-field 
value 6 ~ I \2d\ and than the directed-percolation value 
(0.159). Possibly novel critical behaviour of our model 
is confirmed with calculation of other critical exponents, 
but detailed presentation of our results will be given else- 
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We also examined the d = 2 version of our model. 
Obtained results [l^ suggest that the values of critical 
exponents are very close to those in the d — 1 version. 
The d = 2 version of our model is certainly more realistic 
in the context of tumor-growth modeling. We examined 
their growth starting from a configuration containing a 
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FIG. 4. Time dependence of the density of tumor cells Xt{t) 
(log-log scale) for d = 1, A = 0.1 and (from top) p — 0.356, 
0.3545, 0.354, 0.3538, 0.3537, 0.3536, 0.3535, and 0.353. The 
dotted straight line has the slope corresponding to S = 0.6. 
Inset shows the steady-state density of tumor cells Xt as a 
function of p. In the vicinity of the critical point, xt seems to 
decay linearly (/3 = 1). 
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FIG. 5. Different morphologies of d = 2 tumors growing from 
a single tumor cell surrounded by empty lattice sites. Calcu- 
lations were made on 500x500 lattices and p was chosen such 
that the steady-state density in all three cases was approxi- 
mately equal {xt « 0.05). 



single tumor cell and all other sites empty. Simulations 
were performed for various p and A but such that the 
steady-state density of tumor cells Xt was approximately 
the same and equal to 0.05. Numerical simulations show 
that A strongly affects the shape and dynamics of grow- 
ing tumors (FiglS]). While for small A their shape is 
highly irregular, for large A nearly circular shape with 
a well defined boundary is seen. Apparently, in the for- 
mer case the growth is very much affected by stochastic 
fluctuations while in the latter their role is diminished. 

Different morphologies of growing tumors suggest that 
their other characteristics will also strongly depend on 
A. We measured the tumor survival probability starting 
from the same configurations containing a single tumor 
cell and monitored whether tumor cells survived until a 
given (large) simulation time. The simulations were per- 
formed for several values of p and A and the survival 
probability Psurv is plotted as a function of the steady- 
state density of tumors Xt-, that was estimated indepen- 
dently with standard steady-state simulations (Fig. [5]). 
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FIG. 6. Tumor survival probability Psurv as a function of 
steady-state density xt- Initial configuration contained a sin- 
gle tumor cell surrounded by empty lattice sites. Calculations 
were made on 200x200 lattice and evolution of tumors was 
monitored until t — 10*. Calculations were made for sev- 
eral values of p and the corresponding steady-state density of 
tumors Xt was estimated from steady-state simulations, sim- 
ilarly to the d = 1 case shown in Fig. [T] 

Let us notice that some general arguments developed for 
models with absorbing states suggest that in the active 
phase, both Pgurv and the steady-state order parame- 
ter should scale with the same critical exponent /? [20|. 
Thus, one expects an approximately linear dependence 
^ Xt for sufficiently small xt- Numerical results 
confirm that Pgurv increases linearly with xt but they 
also show a strong dependence on A. Indeed, for the 
same Xt, the tumor survival probability for A = 3 and 
10 is much larger than in the A = case. This is yet an- 
other indication of the importance of angiogenic effects 
on the tumor growth. 

In conclusion, we examined a lattice model of tumor 
growth where the positive correlation of nutrient supply 
with tumor cells distribution mimics angiogenic factors. 
Obtained results show that such a correlation shifts the 
location of a tumor extinction and makes their starvation 
to death more difficult. Moreover, it changes the critical 
behaviour of the model, even though as a model with a 
single absorbing state it should belong to the directed 
percolation universality class. Surprisingly, even in the 
d = 1 version the exponent /3 seems to take the classi- 
cal, mean-field value /3 = 1. However, the exponent S 
describing the time decay of the order parameter at crit- 
icality takes a non-classical value 0.60(5) and most likely 
the model represent a novel critical behaviour. Field- 
theory methods were applied to various models with ab- 
sorbing states J2^ and one can hope that also in this case 
they could provide valuable insight. In the d = 2 ver- 
sion such a correlation affects the morphology of growing 
tumors and substantially raises their survival probabil- 
ity. From the statistical mechanics perspective it would 
be desirable to examine the behaviour of our model for 
— 1 < A < as well as in the limit A ^> oo, where the 
model simplifies and nutrient supply is forbidden to occur 
on empty sites. 
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